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We compare the steady state velocity distributions from our three-dimensional inelastic hard 
sphere molecular dynamics simulation for homogeneously heated granular media, with the predic- 
tions of a mean field- type Enskog-Boltzmann equation for inelastic hard spheres [van Noije & Ernst, 
Gran. Matt. 1, 57 (1998)]. Although we find qualitative agreement for all values of density and 
inelasticity, the quantitative disagreement approaches ~ 40% at high inelasticity or density. By 
contrast the predictions of the pseudo-Maxwell molecule model [Carrillo, Cercignani & Gamba, 
Phys. Rev. E, 62, 7700 (2000)] are both qualitatively and quantitatively different from those of our 
simulation. We also measure short-range and long-range velocity correlations exhibiting non-zero 
correlations at contact before the collision, and being consistent with a slow algebraic decay over a 
decade in the unit of the diameter of the particle, proportional to r~''^^"\ where 0.2 < a < 0.3. 
The existence of these correlations imply the failure of the molecular chaos assumption and the 
mean field approximation, which is responsible for the quantitative disagreement of the inelastic 
hard sphere kinetic theory. 



I. INTRODUCTION 

Granular materials are collections of noncohesive 
macroscopic dissipative particles and are encountered in 
nature and in the industry Q. These materials exhibit 
a wide variety of phenomena depending on the exter- 
nal forcing. The rapid granular flow regime, where the 
collisions are modeled as instantaneous binary inelastic 
collisions, is reminiscent of a gas of hard spheres. Thus, a 
common theoretical approach for this regime as a first or- 
der approximation is to model the system by means of the 
kinetic and the continuum equations for smooth inelas- 
tic hard spheres with a velocity-independent coefRcient 
of restitution In the kinetic theory approach, the 
mean field-type Boltzmann or Enskog-Boltzmann equa- 
tion for inelastic hard spheres is used, and most tech- 
niques are directly transposed from the kinetic theory of 
a gas of elastic hard spheres Q. It is known that this 
formulation is a reasonable description for nearly elastic 
(1 — ^ 1, where e is the normal coefficient of restitu- 
tion) and dilute cases. However, these theoretical models 
include approximations such as the truncation of series 
expansion or hierarchy, and the introduction of equation 
closure. The dissipative nature of the collision modifies 
the physics in a nontrivial way, and the accuracy and the 
limitation of the mean field-type kinetic description with 
these approximations is not yet known. The extension of 
the theory for the more inelastic and dense case, includ- 
ing surface friction, is one of major goals of the current 
inelastic kinetic theory Q] . 

There is an attractor in the phase space of the gran- 
ular media, because inelastic collisions dissipate kinetic 
energy; in the absence of an external energy source, a 
granular medium loses its kinetic energy through colli- 
sions and becomes a static pile. To reach a steady state or 
an oscillatory state, a system of granular media requires 



an external energy source. In this paper, we investigate 
the steady state velocity distributions and velocity corre- 
lations of spatially homogeneous granular media subject 
to a volumetric Gaussian white noise forcing. This sys- 
tem has been studied by several authors |||,^,|l2|,^ as a 
reference system for the kinetic theory of granular media. 
In this system, particles collide inelastically, and execute 
Brownian motion between collisions; the motion is anal- 
ogous to Brownian dynamics of hard sphere suspension, 
but here the kinetic energy is dissipated due to the in- 
elastic collisions rather than hydrodynamic drag. This 
system is far from equilibrium, and the steady state ve- 
locity distribution deviates from the Maxwell-Boltzmann 
(MB) distribution. The velocity distribution of this sys- 
tem was first theoretically studied by van Noije tt al. 
who obtained the approximate solutions to the inelastic 
hard sphere Enskog-Boltzmann equation with a Gaus- 
sian white noise forcing. Their results were tested against 
the Direct Simulation Monte Carlo (DSMC) method of 
the inelastic Enskog-Boltzmann equation j|] , and a good 
agreement was found; it confirmed the accuracy of the 
approximate analysis of van Noije et al., because the va- 
lidity of the DSMC method relies on the validity of the 
inelastic (Enskog-)Boltzniann equation Q]. 

In the current study, we use a large molecular dynamics 
(MD) simulation (using up to 477, 500 particles) of inelas- 
tic hard spheres to investigate the steady state velocity 
distributions, and quantitatively examine the accuracy of 
the inelastic Enskog-Boltzmann equation model for this 
system. Our method is free from the assumptions used in 
the inelastic kinetic theory, and it has an advantage that 
correlations can develop but has a disadvantage that it 
may have a finite simulation size effect. 

We also compare our results with theoretical predic- 
tions of Carrillo et al. ^, who obtained the steady 
state velocity distribution by using the pseudo-Maxwell 
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molecule model 0. The Maxwell molecule model is 
the special case of the inverse power law model for the 
inter-particle potential, whose potential has the form of 
V{r) ~ where r is the inter-particle distance. This 
model has been widely used for analytical studies as a 
reference system for more realistic systems, because the 
model facilitates calculations involving linearized Boltz- 
mann operator. The pseudo-Maxwell molecule model is 
an inelastic analog of the Maxwell molecule model. 

The remainder of the paper is organized as follows. 
Section II presents the method of the numerical sim- 
ulation, and Section III briefly reviews the theoretical 
predictions of the inelastic hard sphere model and the 
pseudo-Maxwell molecule model. In Section IV, simu- 
lation results of the velocity distributions are presented, 
and these are compared with the theoretical predictions 
(Sec. IV A). The simulation results of the velocity cor- 
relations are also presented in this Section (Sec. IV B). 
The paper is summarized and concluded in Section V. 

II. NUMERICAL SIMULATION 

We simulate an ensemble of inelastic hard sphere parti- 
cles of diameter cr in a three-dimensional (3d) cubic box 
of each side 105. 3(t, which are subject to a volumetric 
Gaussian white noise forcing. Particles obtain the kinetic 
energy from the white noise forcing, execute Brownian 
motion in between collisions, and dissipate the kinetic 
energy through inelastic collisions. To be consistent with 
the theoretical studies presented in Section III, we im- 
plement a velocity-independent coefficient of restitution 
e, and neglect the rotational degrees of freedom. We use 
an event-driven MD code originally developed to simu- 
late the patterns in vertically oscillated granular layers. 
Excellent agreement was found between simulations and 
experiments pO| ] . Wc modify this code to implement the 
Gaussian white noise forcing. Brownian dynamics for 
hard sphere simulation was originated in the work of Er- 
mak et al. pT[ |, and the granular analog was introduced 
by Williams et al. for their study of a one-dimensional 
system. 

To implement Gaussian white noise forcing in our MD 
code, we start from a stochastic equation of motion for a 
particle. The equation of motion is given by 

xdt) = :Ft\t) + r,{tl (1) 

where the mass of the particle is the unity, Xi is the i*^ 
cartesian component of the position, T^'^^ is the i*'* com- 
ponent of the forcing due to the collisions, and ri{t) is the 
i*^ component of the stochastic forcing. The stochastic 
forcing term satisfies 

< r,(t) >= 0, (2) 

which assures that the fluctuation cancels out on the av- 
erage, where <> is an ensemble average, and 



< T,{t)T,{t') >= 2F5,j5{t - t'), (3) 

where F is the strength of the correlation, Sij is the Kro- 
necker delta, and 6{t) is the delta function. Eq.(3) assures 
that the collisions well separated in time are statistically 
independent. We assume that all higher-order moments 
of the random variable Ti(t) can be expressed in terms 
of the second moment, which is identical to the assump- 
tion that ri{t) is distributed according to the Gaussian 
distribution [ p^ . We implement an equation equivalent 
to Eq.(l) for each particle with a discrete time interval of 
a fixed size. It can be shown that the discrete Langevin 
equation for the velocity in between collisions subject to 
a Gaussian white noise is given by 

v,{t + dt) =v,{t) + \/Fy/6iG{0,l), (4) 

where F is the same quantity as in Eq.(3), 6t is the time 
interval between the white noise forcing, and G(0, 1) is 
the unit Gaussian distribution random variable. 

In the simulation, Nk particles are randomly chosen 
at every 6t and are kicked in accordance with Eq.(4). 
To avoid the development of a mean flow in the system, 
particles to be kicked are randomly picked pairwise, and 
particles of this pair are kicked in opposite directions with 
the same speed to conserve the momentum. In an event- 
driven simulation, the random kickings are computation- 
ally expensive discrete events, because the collision list 
needs to be updated at each kicking. For the efficiency 
of the simulation, the mean kicking frequency needs to 
be minimized, or the mean kicking time for each particle 
St {= 5t ■ N/Nk) needs to be maximized, where N is the 
total number of particles. We checked in the simulation 
that the velocity distributions do not change as far as St is 
less than l/{5h'coii), where Vcoii is the mean collision fre- 
quency. We also checked that the velocity distributions 
do not depend on the choice of Nk and St in these cases. 
Thus we keep St ^ l/{Wi^coii) throughout the simula- 
tions. The volume fraction varies from i^o{= 4.29%) to 
5j/o(= 21.4%), which corresponds to 95,495 to 477,500 
in the number of particles. Since the collisions are in- 
stantaneous, there is only one time scale determined by 
the granular temperature, ct/a/T, where T is defined in 
Eq.(6). Therefore the temperature only rescales the time, 
and the results are independent of T. However, we fix the 
granular temperature at approximately the same value 
throughout the simulations. We prepare the initial con- 
ditions by locating particles in a regular lattice. We heat 
them, and wait until the transients decay away, ensuring 
that a steady state is reached. The data are taken pe- 
riodically with a fixed time interval At. To assure that 
each data set is statistically uncorrelated. At is chosen 
larger than ^/vcoii in all cases, and At ~ lO/vcoii in most 
cases. We obtain 50 such data sets for each simulation, 
and the error bars appearing in this paper are the stan- 
dard deviations of these data sets. A periodic boundary 
condition is imposed in all direction. 
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III. REVIEW OF THEORY 



The first three Sonine polynomials are given by 



In this Section, we briefly review the results of the 
previous theoretical studies of van Noije et al. [|| and of 
Carrillo et al. |^. Both studies are based on the mean 
field-type inelastic Enskog-Boltzmann equation. 

Following standard procedures of the kinetic theory, it 
can be shown that the inelastic Enskog-Boltzmann equa- 
tion for a system of spatially homogeneous granular par- 
ticles subject to a white noise forcing is given by 



f = Q(/,/) 



Lppf, 



(5) 



where / = /(c, t) is the single particle distribution func- 
tion, Q{f, /) is the collision operator, c — v/ ^j2T{t) 
is the velocity scaled by the characteristic velocity, and 
the mass is the unity. T{t) is the granular temperature, 
defined as the variance of the velocity distribution per 
degree of freedom. For a 3d system. 



< |v(t)- < v(i) > p > 



(6) 



For the inelastic hard sphere model, van Noije et al. 
obtained the steady state solutions for 2d and 3d cases, 
by using the moment method |]l5| which was used in the 
analysis of the velocity distributions of homogeneously 
cooling of a freely evolving system. They expanded the 
solution in Sonine polynomials and neglected terms of 
higher order than the first nonvanishing correction to the 
MB distribution. The steady state solution for a 3d case, 
f^{c), is obtained as 



(10) 



where 



The Fokker-Planck operator Lpp is a diffusion operator 
in the velocity space, which is given by 



L 



FP 



(7) 



where F is the strength of the correlation of Gaussian 
white noise stochastic forcing (defined in Eq.(3)), and 
is the Laplacian in the velocity space. The collision op- 
erator, Q{f, /), is chosen depending on the inter-particle 
potential model, and the most obvious one for the granu- 
lar particles is the inelastic hard sphere collision operator. 

Van Noije et al. ||] obtained the steady state solutions 
to Eq.(5) using the inelastic hard sphere collision opera- 
tor, and Carrillo et al. |8| obtained the steady state so- 
lutions using the pseudo-Maxwell collision operator. In 
both analyses, the Sonine polynomial expansion method 
was used to construct the solution. Sonine polynomials 
are associated Laguerre polynomials and have been used 
to construct solutions to the Boltzmann equation since 
Burnett [ p^ introduced them in the study of nonuniform 
gases. They are exact eigenf unctions of the linearized 
Boltzmann equation for Maxwell molecules. This ex- 
pansion also leads to rapidly converging solutions of the 
linearized Boltzmann equation for other short-range re- 
pulsive potentials The Sonine polynomials of lower 
parameter 1/2 is defined by 



(i+n)! 



i-xr. 



(8) 



In the current study, is the argument of Sonine poly- 
nomials. The orthogonality relation with this argument 
is 



(9) 



Imb (c) 



=c^ exp(-c^) 



(11) 



is the MB distribution (multiplied by an integration fac- 
tor), and 



16(l-e2)(l-2e2) 



73 + 56rf - 24ed - 105e + 30(1 - e)e^ ' 



(12) 



where d is the dimensionality of the system (d = 3 in the 
current study), and HS stands for inelastic hard spheres. 

For the pseudo-Maxwell molecule model, Carrillo et 
al. Q obtained the steady state solution to Eq.(5) by 
doubly expanding the solution in the energy dissipation 
rate, e = (1 — e^)/4, and in Sonine polynomials. The de- 
viation from the MB distribution was obtained up to the 
order of . The solution with the first nonvanishing 

correction to the MB distribution is of the order of e^, 
which is given by 



/(MAf)(c) = /mb(c)[1 



(13) 



where MM stands for pseudo-Maxwell molecules. The 
first order deviations from the MB distribution in the 
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two models have the same basis function s[y^{c^), 
second Sonine polynomial. We rewrite the normalized 
deviations from the MB in the two models as follows, 



Af"^ie,c) 

/mb(c) 
A/^^(e,c) 

fMB{c) 



(2) 



arHe)S[%{c% 



(14) 
(15) 



where 
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af^(e) = (l- 6^)74 = 4e2. (16) 

The coefficients of the second Sonine polynomial for the 
two models, a2^{e) and a^^{e), are compared in Fig. 
1. 

To summarize, the predictions of the inelastic hard 
sphere model and those of the pseudo-Maxwell molecule 
model are different mainly in two ways: 

(1) There is a crossover from positive to negative 
values in a|^'^(e) as e increases, and it is negative for 
l/^/2 < e < 1.0, while af^(e) is positive definite, 
4e^. The high energy tail is always overpopulated in the 
pseudo-Maxwell molecule model, while it is underpopu- 
lated for e > 1/a/2 in the inelastic hard sphere model, 
when the series is truncated at this order. 

(2) The deviation from the MB distribution in the 
pseudo-Maxwell molecule model is quantitatively much 
larger than that in the inelastic hard sphere model. Note 
that the analysis of the pseudo-Maxwell molecule model 
uses a small inelasticity assumption, while the analysis of 
the inelastic hard sphere model has no such assumption. 

IV. SIMULATION RESULTS 
A. velocity distributions 

We measure the velocity of each particle periodically in 
time with a fixed time interval of At ^ ^Q/vcoiu which is 
chosen to assure that each data set is statistically uncor- 
related. The measured velocities are binned, and the bin 
size is Ac = 0.1. The velocity distributions for two coeffi- 
cients of restitution, e = 0.1 and e = 0.9, are obtained in 
the simulation (Fig. 2). As in the inelastic hard sphere 
theory, high velocity tail is overpopulated for e = 0.1, 
and slightly underpopulated for e = 0.9, compared to 
the MB distribution. 

We define the normalized deviation from the MB dis- 
tribution, g{c), which is obtained in the simulation, 

/(md)(c)=/mb(c)[1+5(c)], (17) 

where the subscript (MD) means this is obtained in the 
MD simulation. Assuming g{c) is a smooth function in 
the scale of Ac (<C 1), g{c) can be approximated as 

g(^o + ^)^ (18) 

J Co -/tT 

The numerator in Eq.(18) is the total number of particles 
in the bin, which is a number counted in the simulation, 
and the denominator can be evaluated using the error 
function table. g{c) calculated using Eq.(18) for the data 
in Fig. 2 is shown in Fig. 3. For direct comparison 
between our simulation results and the above theoreti- 
cal predictions, we calculate the coefficient of the second 
Sonine polynomial, a^^, from our measurements. As in 



the theoretical analyses, we assume f^j^jj-^ (c) is expanded 
in Sonine polynomials, 

oo 

flMD){c)=fMB{c)[l + Y,<''S'^,{c^)], (19) 
fe=2 

where o/^^^ is not included, because it is identically zero 
in theory. We checked in the simulation that af^^ is less 
than 10~^ in all cases. 

In the simulation, we can use either of the following 
two formulae to calculate , the coefficient of the k*^ 
Sonine polynomial. Firstly, a^^'s can be obtained from 
the following numerical integration, using the orthogo- 
nality relation for Sonine polynomials, 

«r = f flMD) (c)5i5^(c^)rfc. (20) 

Secondly, we can make use of the recurrence relation of 
a^^'s. Starting from the definition of Sonine polynomi- 
als, Eq.(8), it can be shown that for fc > 2, satisfies 

-(2fc + l)!!<'^ >+i V ' W 

(21) 

where < c^*^ > is the 2k*^ moment. It is straightforward 
to numerically evaluate the integration in both formulae. 
The above two relations, Eq.(20) and Eq.(21), are math- 
ematically identical. However, the errors of the lower k 
coefficients are accumulated on higher k coefficients in 
Eq.(21), while each is determined independently in 

Eq.(20). The results obtained using Eq.(20) and Eq.(21) 
are nearly the same for small k, and they deviate more 
for larger k. of^^ is obtained using Eq.(20), which are 
compared with the theoretical predictions Eq.(12), a2^, 
and Eq.(16), ^ (Fig. 4). The simulation results devi- 
ate more from the predictions of the inelastic hard sphere 
kinetic theory as the system becomes more inelastic, and 
has a crossover between the positive and the neg- 
ative values at around e ~ 0.8, while it was predicted 
to occur at l/\/2 by the inelastic hard sphere theory. 
The pseudo-Maxwell molecule model does not predict a 
crossover behavior of al^*'^, and the predictions of this 
model deviate quantitatively from the kinetic theory and 
simulations of the inelastic hard sphere model. 

In the simulation, all aj^'^'s can be obtained by using 
Eq.(20). For larger k, the weight of the high velocity data 
exhibiting strong fluctuations become larger, and so does 
the uncertainty of . We calculate up to for var- 
ious coefficients of restitution (Fig. 5). For nearly elastic 
case, for instance e = 0.9, the series converges rapidly 
as in many cases of elastic hard spheres. However, as 
the system becomes more inelastic, the series converges 
more slowly. For e = 0.1, af ^ is about 30% of a^^. It 
is shown that how the first few nonvanishing coefficients 
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of the Sonine polynomial series obtained in the simula- 
tion affect the fitting (Fig. 6). For e = 0.1, a^^^ fits 
the measured data well up to w ^ 3\/2r ^ 4a/T, but 
the high energy tail is fitted well only when the series 
are kept up to a*^^. For e = 0.9, it shows the same 
tendency. When we fit the data only with in this 

case, /(\/£))(c) becomes negative for c > 4.1 and becomes 
negative infinity as c goes to infinity. Such nonphysical 
behavior disappears when we include a^^^ . 

In the theoretical predictions, the steady state velocity 
distribution and a|^'^ do not depend on the density, and 
it is a function of the coefficient of restitution and the di- 
mensionality only. However, we find that it also depends 
on the density (Fig. 7). The value of a2^^ decreases 
as the system becomes more dilute. For the smallest 
volume fraction we used, i/q, the value of a2^^ deviates 
by only few percent from the predictions of the inelastic 
hard sphere theory. We plot a|^^'s for various densities 
in Fig. 8 (af^^'s are very small, as shown in Fig. 5). It 
follows the same tendency as Oj^^ does; as the system 
becomes denser, af^^ gets larger. 

Finally, we examine the asymptotic behavior of the 
high energy tails of the velocity distribution functions. 
Van Noije et al. ||] found for this system a high energy tail 
shows an asymptotic behavior ~ exp(— ^'c^^^). We mea- 
sure the velocity distribution function f'^MD)('^)^ which is 
defined as 

1 = J /(MD)(c)c?c 
= y"/(MD)(c)4^c2dc 
= / fiMD)ic)dc. (22) 

To investigate the power of the argument of the expo- 
nential function, we renormalize /(\f (c) with its max- 
imum value. We observe the crossover behavior from 
~ exp(— ^c^) to ~ exp(— y^'c"^/^) as c increases, for 
e < 0.5 (Fig. 9). 

B. Velocity correlations 

Two major approximations imposed in the kinetic the- 
ory discussed in Section III are the mean field approxima- 
tion and the truncation of the hierarchy by introducing 
the molecular chaos assumption. In this Section, we ex- 
amine the validity of each of the above approximations 
for this system by quantitatively investigating the par- 
allel velocity correlations, in long-range and short-range 
respectively. It is known that a system of granular me- 
dia exhibits strong spatial correlations |jl^,|l^,|l8| , such 
as velocity correlations, but there is still the lack of a 
quantitative study of a 3d system. We suggest that the 
deviation of Oj^^ from originates from the failure of 
the above approximations. 

We define the parallel velocity correlation function as 



< Ci,||C2,|| >= ;^ J < c(r + r') • fp(r + r')c(r') • fp(r') > dr' , 

(23) 

where f = r/|r|, and p(r) is the local particle density, 

N 

p(r)=^^(r-r,). (24) 

i=l 

We approximately evaluate this quantity using the fol- 
lowing formula, 

<^i.lF2^ll>-E^l^' (25) 

where the parallel direction is along the line of centers of 
the particle pair under consideration, Nr is the number of 
particles in a shell of thickness Sr and inner radius r, and 
the summation is done over Nr. We use Sr = 0.10530-. 



1. Long-range correlations 

The parallel velocity correlations are obtained by av- 
eraging over 50 statistically uncorrelated data sets. We 
calculate those for various coefficients of restitution as a 
function of dimensionless distance r/a (Fig. 10), and for 
various densities (Fig. 11). The data are shown only for 
r/<T < 30, because they are subject to the finite system 
size effect for larger r/a. The parallel velocity correla- 
tions in our simulation are consistent with slow algebraic 
decay over a decade, ~ 7-~(i+'5)^ where 0.2 < S < 0.3. 
This behavior is close to the theoretical prediction of van 
Noije et al. |l^, who predicted the r~^ power law from 
the mode coupling theory. The correlations in our simu- 
lation get stronger for more inelastic or more dense sys- 
tem, which implies that the mean field approximation is 
reliable only for nearly elastic and dilute cases. 



8. Short-range correlations 

In this section, we investigate the velocity correlations 
at contact before the collision, to examine the validity 
of the molecular chaos assumption. A non-zero value 
of these correlations is the signature of the failure of 
the molecular chaos assumption; the velocities are more 
"parallelized" after the collision, since only the normal 
component of relative velocity of colliding particles are 
reduced at the collision, which means that if the veloci- 
ties are already parallelized before the collision, it would 
indicate that the colliding pair have "memory" on the 
collisions in the past, and they are correlated. We find 
that for high inelasticity and density, the pre-coUisional 
parallel velocity correlation value reaches up to ~ 15% of 
the temperature. Even for dilute or nearly elastic cases, 
these are not negligibly small. 
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We calculate the velocity correlations of pre- and post- 
coUisional states by evaluating Eq.(25) for approaching 
(ri2 • Ci2 < 0) and separating particles (ri2 • C12 > 0) re- 
spectively, where X12 = Xi — X2 for x = r or c. We calcu- 
late the pre-collisional state of short-range (1 < r/cr < 2) 
parallel velocity correlations (Fig. 12), and their post- 
collisional state (Fig. 13), for various coefficients of resti- 
tution. The values of the pre-collisional correlations are 
not negligible compared to the temperature of the sys- 
tem. Post-collisional correlations are more than twice as 
large as pre-collisional correlations. 

The maximum values of the velocity correlations in 
Fig. 12 and in Fig. 13 are not the values at the con- 
tact, r = a, because of the finite size of the bins in the 
measurements; instead, those are values at r/a = 1.053. 
The values at r — a are estimated by extrapolating the 
data in the interval 1 < r/cr < 2 using the least square 
fit with fifth-order polynomials. We estimate the velocity 
correlations at contact, r = cr, as a function of the coef- 
ficient of restitution (Fig. 14), and as a function of the 
density (Fig. 15). Since the velocity correlation varies 
rapidly as r/a decreases to 1, these estimations may be 
regarded only as approximate lower bounds. The velocity 
correlations at contact in this system show almost linear 
behavior both in density and the coefficient of restitution. 



V. DISCUSSION 

We have investigated the velocity distributions and 
parallel velocity correlations of 3d homogeneously heated 
granular media for various densities and inelasticities, us- 
ing an inelastic hard sphere MD simulation. The devia- 
tions from the MB distribution in our simulations quali- 
tatively agree with the results of the mean field-type in- 
elastic hard sphere kinetic theory but we found that 
there is systematic quantitative difference. 

We observed the high energy tails are consistent with 
^ exp(— yi'c^/^) for e < 0.5, but not for others. Since the 
elastic case (e — 1.0) has no crossover from ~ exp(— y^c^) 
to ^ exp(— ^'c^/^), we expect that this crossover behav- 
ior may occur at higher velocities as e approaches to 1.0, 
if it occurs. However, we were not able to check whether 
the crossover occurs for e > 0.5 or never occurs, because 
our system is finite. It is interesting to note that the 
same behavior was experimentally observed in a system 
with different forcing mechanism ]l9t . 

We found that the steady state velocity distributions 
in the simulation depend on the density as well as the 
coefficient of restitution, while they depend only on the 
latter in the theory. The discrepancy between our sim- 
ulation results and the theoretical predictions increases 
as the system becomes more inelastic or more dense, and 
the quantitative disagreement reaches up to ~ 40%. This 
behavior is consistent with the results of van Noije et 
at |l^,|2^, who found that the collision frequency mea- 
sured in a 2d MD simulation deviates more from the pre- 



dictions of the inelastic Enskog-Boltzmann equation as 
the system becomes more inelastic. We suggest that the 
disagreement originates from the failure of two major ap- 
proximations in the theory, the mean field approximation 
and the molecular chaos assumption. To examine the ac- 
curacy and the limitation of these approximations, we 
quantitatively investigated the parallel velocity correla- 
tions of this system. 

We found that the long-range parallel velocity cor- 
relations are consistent with a slow algebraic decay, ~ 
where 0.2 < 6 < 0.3. This result is close to the 
theoretical predictions of van Noije et al. ||l^, who renor- 
malized various quantities such as the collision frequency 
using the mode coupling theory and predicted power 
law for velocity correlations in the system we studied. 
Because of these strong correlations, the mean field ap- 
proximation is not a good one unless the system is nearly 
elastic or very dilute. 

We also found that the velocity correlations at con- 
tact before the collision are not negligible. We measured 
the short-range velocity correlations of pre- and post- 
collisional states separately to examine the validity of the 
molecular chaos assumption. The pre-collisional correla- 
tions at contact are about a half of the post-collisional 
correlations. The correlations at contact are almost lin- 
early proportional to both the density and the coefficient 
of restitution, which is consistent with the recent results 
of Soto et al. who studied the velocity correlations 
of a 2d homogeneously cooling granular media in nearly 
elastic regime. 

We also examined the convergence of the Sonine poly- 
nomial expansion technique used in the inelastic kinetic 
theory, and found that the series converges more slowly 
as the system becomes more inelastic. 

Finally, we compared the steady state velocity distri- 
butions in the simulations with the theoretical predic- 
tions of Carrillo et al. , who studied the current system 
using the pseudo-Maxwell molecule model. We found 
that the velocity distribution function predicted by this 
model differ qualitatively from those predicted by the in- 
elastic hard sphere model. 
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A. Derivation of Eq. (18) 

The deviation from the MB distribution, g{c), is de- 
fined as in Eq.(17), 



fiMD){c) = fMB{c)[l+g{c)]. 



(Al) 



In the simulation, the number of particles in each bin is 
measured; 



/■Co+Ac /■ 
/ f(MD){c)dc= / 



Co+Ac 



(MD){c)dc= I -^e " c [l + g{c)]dc, 



(A2) 

where the bin size Ac is assumed to be very small. It 
is possible to numerically solve for g{c) from Eq.(A2), 
however, we get an approximate expression for g{c) by 
assuming g{c) is smooth in the seal of Ac; 



Ac. r-.+^c 



J e-" c'g{c)dc « g{co + ^) J 



e ° c^dc. 



Eq.(A2) can be read 



9{co + -tt) 



(MDy 



2' - po+Ac 4 g_e^^2^^ 
J Co v'"' 



(A3) 



- 1- (A4) 



B. Derivation of Eq. (21) 



Starting from the definition of Sonine polynomials, 



!L (1+n)! 
'-r'a{k+Py-in-p)\p\ 



{-xf, (Bl) 



p=o '^2 
it can be shown that 

_ ' (^1)^+^(1 + fc)!fc! (,_,) 

We assume the following Sonine polynomial expansion of 
the distribution function, 



= /mb(c)[1 + £ ^5^^) (c^)]. (B3) 
Using Eq.(B2) and Eq.(B3), 2k*^ moment reads 



(?^[i+x:(-i)'=+^Qan], (B4) 

p=0 



_(-lf2_ 
(2A; + 1)!! 



< c^" > +(-1) 



fe+i 



MD 
k—p ' 



(B5) 



where k > 2. 



or 
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FIG. 1. The deviation from the MB distribution predicted by the theory, 02, as a function of the coefficient of restitution e. 
a2^{e) is the coefHcient of the second Sonine polynomial from the inelastic hard sphere model, and 02"^' (e) is that from the 
pseudo-Maxwell molecule model. 
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FIG. 2. The steady state velocity distributions, f(MD)ic), obtained from the simulation for two coefficients of restitution, 
e = 0.1 and e = 0.9. The volume fraction is 5vo- The thick solid line is the MB distribution function. 
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FIG. 3. The normalized deviation from the MB distribution, g{c), obtained from the simulation for two coefficients of 
restitution, e = 0.1 and e = 0.9. The volume fraction is Sz/q. The error is large for c > 3, because it involves division by a very 
small number. 
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FIG. 4. A comparison of a2's from the three models, the kinetic theory of the inelastic hard sphere model (solid curve), 
the kinetic theory of the pseudo-Maxwell molecule mode (dashed curve), and the MD simulation (open circles). The volume 
fraction is 5vo for the simulation. 
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FIG. 6. The normalized deviations from the MB distribution, g{c), obtained in the simulation (open circles), are compared 

with the fitted curves using the coefficients of the Sonine polynomial expansion calculated using the simulation data, which 
are shown for two coefficients of restitution, e = 0.1, and e = 0.9. The terms are successively added to the Sonine polynomial 
expansion up to a^^. The volume fraction is for both, and the error bars for the data are not included for better comparison 
with fitting curves. 
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FIG. 7. The values of (12'^ obtained in the simulation as a function of density for two coefficients of restitution, e — 0.1 
and e = 0.5. Dashed lines are the predictions of the inelastic hard sphere kinetic theory, which do not depend on the density. 
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FIG. 8. The density dependence of for two coefficients of restitution, e = 0.1 and e = 0.5, obtained in the simulation. 

Predictions arc not available, because they have not been calculated in the theoretical studies. Error bars for these data are 
bigger than those for af^^ (Fig- 7), because in calculating a^'^ , higher velocity data have more weight exhibiting stronger 
fluctuations. 
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FIG. 9. A crossover behavior of the velocity distribution function f^MD){'^) (open circles) from exp(— ^c^) (compared 
with dashed lines) to ~ exp{—A'c^^^) (compared with solid lines) is observed for higher inelasticities (e < 0.5). The volume 
fraction is buo, and A and A' are arbitrary constants. 
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FIG. 10. The parallel velocity correlations for various coefficients of restitution. The volume fraction is 5fo. Dashed lines 
are the curves following the power law, which are included for comparison. The curves deviate from the power law for r/a> 20, 
because of the finite system size effect. Error bars are not shown for clarity. 
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FIG. 12. The short-range pre-colhsional parallel velocity correlations for various coefficients of restitution. The volume 
fraction is 5z/o. 
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FIG. 14. The pre- and post-collisional parallel velocity correlations at r/cr = 1.053 (solid lines) and estimated values at 
contact, r ja = 1 (dashed lines), as a function of the coefBcient of restitution. The volume fraction is Bfo- The values at contact 
are obtained by extrapolating the data in Fig. 12 and Fig. 13, using fifth order polynomials. These are extrapolated values 
from the average values, and the error bars are not systematically determined. The error may be of the same order as the 
values at r/a = 1.053 in Fig. 12 and Fig. 13. 
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FIG. 15. The pre- and post-coUisionaJ parallel velocity correlations at r/a = 1.053 (solid lines) and estimated values at 
contact, r/a = 1 (dashed lines), as a function of the density. The coefficient of restitution is 0.1. 
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